Associations between DNA methylation and gene regulation depend on chromatin accessibility during transgenerational plasticity

Background Epigenetic processes are proposed to be a mechanism regulating gene expression during phenotypic plasticity. However, environmentally induced changes in DNA methylation exhibit little-to-no association with differential gene expression in metazoans at a transcriptome-wide level. It remains unexplored whether associations between environmentally induced differential methylation and expression are contingent upon other epigenomic processes such as chromatin accessibility. We quantified methylation and gene expression in larvae of the purple sea urchin Strongylocentrotus purpuratus exposed to different ecologically relevant conditions during gametogenesis (maternal conditioning) and modeled changes in gene expression and splicing resulting from maternal conditioning as functions of differential methylation, incorporating covariates for genomic features and chromatin accessibility. We detected significant interactions between differential methylation, chromatin accessibility, and genic feature type associated with differential expression and splicing. Results Differential gene body methylation had significantly stronger effects on expression among genes with poorly accessible transcriptional start sites while baseline transcript abundance influenced the direction of this effect. Transcriptional responses to maternal conditioning were 4–13 × more likely when accounting for interactions between methylation and chromatin accessibility, demonstrating that the relationship between differential methylation and gene regulation is partially explained by chromatin state. Conclusions DNA methylation likely possesses multiple associations with gene regulation during transgenerational plasticity in S. purpuratus and potentially other metazoans, but its effects are dependent on chromatin accessibility and underlying genic features. Supplementary Information The online version contains supplementary material available at 10.1186/s12915-023-01645-8.


Table of Contents:
A. General supplemental figures (p. 3 -8) B. Specification and diagnostics for selected model of differential expression under maternal upwelling as a function of intron differential methylation (p. 9 -15) C. Specification and diagnostics for selected model of differential exon use under maternal upwelling as a function of differential exon methylation ( Figure S6: Differential intron methylation affected expression conditional upon intron length. Differential gene expression under maternal upwelling is plotted against mean intron differential methylation. Data are grouped based on genic intron length quartiles. 'Long introns' and 'short introns' denote highest and lowest quantiles. Linear regressions are fitted across observed values. Average logFC across binned intron differential methylation is plotted as black points ± SE. n = 12 RNA-seq and RRBS replicate libraries.

General supplemental figures
Specification and diagnostics for selected model of differential expression under maternal upwelling as a function of intron differential methylation Figure S7: Posterior distributions and MCMC chains for β parameters of the selected model predicting differential expression under maternal upwelling as a function of differential intron methylation. "_Z" is appended to the end of parameters that were scaled to Z-scores during model fitting in order to improve run time and convergence of MCMC chains. 'GE_logCPM' denotes logCPM of gene expression. 'int_meth' denotes baseline % methylation of intronic CpGs. 'int_meth_logFC' represents log2FC of differential intron methylation across genes. 'TSS_dens_norm' represents the density of chromatin accessibility at ± 500 bp TSS regions. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries.
Specification and diagnostics for selected model of differential expression under maternal upwelling as a function of intron differential methylation Figure S8: Leave-one-out estimates of leverage for observed data fit to selected model of differential expression under maternal upwelling as a function of differential intron methylation. Observations with pareto shape k > 0.4 are deemed to have moderate leverage capable of biasing model fitting. Observations with pareto shape k > 0.7 possess high leverage. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries.
Specification and diagnostics for selected model of differential expression under maternal upwelling as a function of intron differential methylation Figure S9: Correlation matrix of β posterior draws for fixed effects in selected model of differential expression under maternal upwelling as a function of differential intron methylation. Darker blue depicts greater point density. "_Z" is appended to the end of parameters that were scaled to Z-scores during model fitting. 'GE_logCPM' denotes logCPM of gene expression. 'int_meth' denotes baseline % methylation of intronic CpGs. 'int_meth_logFC' represents log2FC of differential intron methylation across genes. 'TSS_dens_norm' depicts the density of chromatin accessibility at ± 500 bp TSS regions. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries.
Specification and diagnostics for selected model of differential expression under maternal upwelling as a function of intron differential methylation Figure S10: Posterior predictive check of selected model predicting differential expression under maternal upwelling as a function of differential intron methylation. The x-axis depicts Z scorescaled differential expression logFC values. The y axis the density distribution of observed and predicted logFC. The black line (γ) depicts the distribution of observed data. Blue lines (γrep) depict iterative distributions of model predictions. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries. Specification and diagnostics for selected model for differential exon use under maternal upwelling as a function of differential exon methylation Figure S13: Posterior distributions and MCMC chains for β parameters of the selected model predicting differential exon use under maternal upwelling as a function of differential exon methylation. "_Z" is appended to the end of parameters that were scaled to Z-scores during model fitting in order to improve run time and convergence of MCMC chains. 'exon_num' represents exon number. 'logCPM' denotes logCPM of gene expression. 'exon_DM_logFC denotes log2FC values of exon differential methylation. 'intron_len' represents total genic intron length. 'exon_dens_norm' depicts the density of chromatin accessibility at across exons of the associated gene. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries.  Table S2: Posterior intervals representing probability of direction tests applied to selected model of differential exon use under maternal upwelling. Significant effects have posterior probabilities for which >95% of the distribution falls above or below 0. "_Z" is appended to the end of parameters that were scaled to Z-scores during model fitting in order to improve run time and convergence of MCMC chains. 'exon_num' represents exon number. 'logCPM' denotes logCPM of gene expression. 'exon_DM_logFC' represents log2FC of differential methylation at single exons. 'exon_dens_norm' depicts the density of chromatin accessibility at all exons within a gene. An asterisk denotes methylation parameters modeled with an error term equaling inverse exon CpG coverage per gene. Asterisks denote significant fixed effects as evidence by probability of direction.

Specification and diagnostics for selected model for differential exon use under maternal upwelling as a function of differential exon methylation
Specification and diagnostics for selected model for differential exon use under maternal upwelling as a function of differential exon methylation Figure S14: Leave-one-out estimates of leverage for observed data fit to selected model of differential exon use under maternal upwelling as a function of differential exon methylation. Observations with pareto shape k > 0.4 are deemed to have moderate leverage capable of biasing model fitting. Observations with pareto shape k > 0.7 possess high leverage. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries.
Specification and diagnostics for selected model for differential exon use under maternal upwelling as a function of differential exon methylation Figure S15: Correlation matrix of β posterior draws for fixed effects in selected model of differential exon use under maternal upwelling as a function of differential exon methylation. Darker blue depicts greater point density. "_Z" is appended to the end of parameters scaled to Zscores during model fitting. 'exon_num' represents exon number. 'logCPM' denotes logCPM of gene expression. 'exon_DM_logFC' represents log2FC of differential exon methylation across genes. 'exon_dens_norm' depicts the density of chromatin accessibility at exons within a gene. "intron_len" represents total intron length of a gene. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries.
Specification and diagnostics for selected model for differential exon use under maternal upwelling as a function of differential exon methylation Figure S16: Posterior predictive check of selected model predicting differential exon use under maternal upwelling as a function of differential exon methylation. The x-axis depicts Z scorescaled differential exon use ΔlogFC values. The y axis the density distribution of obersved and predicted ΔlogFC. The black line (γ) depicts the distribution of observed data. Blue lines (γrep) depict iterative distributions of model predictions. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries.
Specification and diagnostics for selected model for differential exon use under maternal upwelling as a function of differential exon methylation Figure S17: Predictions of differential exon use (DEU) under maternal upwelling by selected model relative to exon differential methylation (DM), exon accessibility (columns), and total genic intron length (rows). Individual points depict fitted values per exon ± 95% credibility intervals. 'Low' and 'high' groupings of exon accessibility and intron length represent observations in the bottom and top quartiles of these variables. Blue solid lines depict fitted regressions to predicted DEU. Red dashed lines depict unfitted regressions to observed DEU. n = 12 RNA-seq and RRBS replicate libraries; n = 3 ATAC-seq replicate libraries.
Specification and diagnostics for selected model for differential exon use under maternal upwelling as a function of differential exon methylation